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The formation mechanism of planetesimals in protoplanetary discs is hotly debated. Currently, the favoured model involves the 
accumulation of meter-sized objects within a turbulent disc, followed by a phase of gravitational instability. At best, one can simulate 
a few million particles numerically as opposed to the several trillion meter-sized particles expected in a real protoplanetary disc. 
Therefore, single particles are often used as super-particles to represent a distribution of many smaller particles. It is assumed that 
small-scale phenomena do not play a role and particle collisions are not modelled. The super-particle approximation is not always 
valid when applied to planetesimal formation because the system can be marginally collisional (of order one collision per particle per 
orbit). The super-particle approximation can only be valid in a collisionless or strongly collisional system, although, in many recent 
numerical simulations this is not the case. 

In this work, we present new results from numerical simulations of planetesimal formation via gravitational instability. A scaled 
system is studied that does not require the use of super-particles. This system is simplified for computational practicality and proper 
identification of important processes: 1) the evolution of particles is studied in a local shearing box; 2) the particle-particle interactions 
such as gravity, physical collisions, and gas drag are solved directly assuming a constant background shear flow without any feedback 
from the particles. We find that the scaled particles can be used to model the initial phases of clumping if the properties of the scaled 
particles are chosen such that all important timescales in the system are equivalent to what is expected in a real protoplanetary disc. 
Constraints are given for the number of particles needed in order to achieve numerical convergence. 

We compare this new method to the standard super-particle approach. We find that the super-particle approach produces unreliable 
results that depend on artifacts such as the gravitational softening in both the requirement for gravitational collapse and the resulting 
clump statistics. Our results show that short-range interactions (collisions) have to be modelled properly. 

Key words, planet formation - planetesimals - W-body simulations - convergence - accretion disc 



1. Overview 

Extrasolar planets have been observed around a variety of par- 
ent stars from pulsars to solar-type stars to M- dwarfs (see e.g., 
IChauvin et ail l2004t IWolszczan & Fraill 1 1992b indicating that 
planet formation is common and successful in a broad range of 
environments. However, the process of planet formation itself 
is not directly observable, leaving theory and numerical simu- 
lations to fill in the blanks between observations of hot circum- 
stellar discs around young stars and planets orbiting mature stars. 
One of the most important unanswered questions in the theory of 
planet formation is what is the mechanism for planetesimal for- 
mation, i.e., the process by which the building blocks of planets 
are formed. 

There are two mai n theories for plane tesimal formation: 
mutual collisions (e.g.. iHavashi et all 1977b and grav itational 
instability in the dust layer ( Goldreich & Ward Il973l) . In the 
first hypothesis, dust particles grow as the result of accretion- 
dominated collisions. Although the formation of planetesimals 
by mutual collisions is consistent with meteoritic evidence, the 
collision speed between dust particles (or aggregates) must be 
much slower than the typical velocity dispersio n in a standard 
proto stellar disc to avoid destructive collisions dBlum & Wurml 
120081) . In addition, the planetesimal formation process is so 
slow that meter-sized particles are in danger of spiraling into 
the star before growing large enough to decouple from the gas 



dWeidenschillind[l977bl) . Even if km-sized planetesimals were 
able to form, they would be in danger of being grou nd down 
again by mutual collisions dStewart & Leinhardl l2009). 

Gavitational instability is often considered to be a solution 
to most of these problems because the intermediate sizes are 
avoided all together. In this theory, the dust layer becomes dense 
enough for the Keplerian shear and velocity dispersion of the 
dust particles to be unable to support the dust against its own 
self gravity. The dust then collapses into clumps that eventu- 
ally cool via drag forces and mutual collisions into planetesi- 
mals. Many different authors have worked on this subject. Until 
recently, the focus has been on quiet, non-t urbulent, and low 
densi t y regions of the accretion disc (see e.g., iMichikoshi et al.l 
I2009L l2007t iTanga et al.ll2004 . However, the turbulent gas in 
the protoplanetary nebula stirs the dust, which increases the ve- 
locity dispersion of the dust particles. Several ideas have been 
proposed to overcome the turbulence-induced mixin g of the dust 
particl es and create localized clumps. For example. ICuzzi et ail 
(120081) suggest that the same turbulence that stirs the dust on 
larger scales may also collect t he dust particles on sm all scales. 
A similar idea was proposed bv lJohansen et alj d2007l) . in which 
dust particles are localized in to clumps using both turbu lence 
and the streaming instability dYoudin & Goodman! 1200 5?). The 
clumps then become gravitationally unstable. A third hypothesis 
suggests that large structures, such as vortices, may be able to 
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collect and protect dust pa rticles from the turbulent background 
dBarge & Sommerialll995l) . In this paper, we focus on the gravi- 
tational collapse in a very dense and turbulent region of the pro- 
toplanetary disc. The overdensity in the dust layer, which is ap- 
proximately two orders of magnitude higher than the standard 
minimum-mass solar nebula, might have been created by any of 
the processes described above. 

Numerical simulations must be used to test models of plan- 
etesimal formation. However, the separation of scales in the 
problem is huge. A meter-sized object is more than 12 orders of 
magnitude smaller than the protoplanetary disc thickness forcing 
numericists to use super-particles. These super-particles have a 
mass much higher than the mass of a single dust particle, but the 
forces (e.g., drag forces ) acting on them ar e equivalent to those 
of a single dust particle. Tanga et ail d2004l) were able to use this 
approximation to succesfully model the gravitational collapse in 
a regime in which collisions are unimportant. However, when the 
surface density is increased by two orders of magnitude, as men- 
tioned above, the system does become marginally collisional and 
therefore the super-particle approach breaks down. 



iMichikoshi et alj d2009l) perform a similar set of simulations, 
also in a low density, non-turbulent region of the disc. Although 
the authors describe their particles as super-particles, the way 
in which they accurately treat the collisions is the same as the 
method presented in this paper. Other approaches use Monte 
Carlo type methods to res olve statistically the outco me of over- 
lapping particles (see e.g. jLithwick & Chiandl2007h . 

In the following, we look carefully at the numerical require- 
ments of modelling gravitational instability accurately and test 
the validity of using super-particles in a high density region. 
For simplicity, we begin with a system without turbulence and 
assume that the surface density has already been increased by 
turbulence or some equivalent process. We find that the super- 
particle approach gives erratic results when collisions become 
important. 

We therefore go on and present a method that does not use 
super-particles. We simulate a scaled system in which the num- 
ber of particles is dramatically lower than in a real protoplane- 
tary disc. To keep the surface density the same, the mass of in- 
dividual particles has to be increased. The drag force produced 
by the surrounding gas is scaled accordingly, so that the stop- 
ping time remains constant. Furthermore, we increase the phys- 
ical size of each particle. This allows us to keep the collision 
timescales in simulations of different particle numbers exactly 
the same. These requirements lead to results that are independent 
of the number of particles, which is the only free parameter in the 
simulation, and we call those simulations converged. The results 
agree with other simulations performed in the context of S aturn's 
rings dGoldreich & Tremaindl982l:lDaisaka^fc Idall999h . In this 
paper, we argue that this method should also be used in current 
simulations of planetesimal formation. 

The outline of this paper is as follows. In Sect. |2l we de- 
fine the important timescales in the problem and show that the 
super-particle approach breaks down in high density regions. In 
Sect. [3] we discuss the numerical method and the initial condi- 
tions of our simulations. Results of numerical simulations using 
both the super-particle approach and our scaling method are pre- 
sented in Sect. [4] We conclude the paper with resolution con- 
straints for numerical simulations and discuss the implications 
for the planetesimal formation process through gravitational in- 
stability in Sect. [5] 



2. Orders of magnitude 

2. 1 . Definition of timescales 

The aim of this work is to simulate the dynamics of dust (or 
particles) interacting gravitationally inside a protoplanetary disc. 
Each particle is subject to five different physical processes, each 
operating on various typical timescales: 

Stopping time T s 

Each particle feels the effects of the surrounding gas 
through a linear drag force. This force can be written as 
Fdrag - "^(Vg ~ v p)> where t s is the stopping time of a 
particle of mass m p , v„ is the velocity of the gas, and v p is 
the velocity of the particle dWeidenschilIing|[l977al) . 

Physical collision timescale T c 

Dust particles will suffer a physical collision with another 
dust particle on a timescale of t c = {cr c v p nf x , where o~ c 
is the geometrical cross-section of the particles, \ p is the 
particle velocity dispersion, and n is the number density of 
particles. 

Orbital timescale T e 

The timescale associated with the particles orbiting the cen- 
tral object in a local shearing patch is the epicyclic period 
T e = 2n£l , where O is the angular velocity at the semi- 
major axis of the shearing box. 

In this paper, we consider two limits for the gravitational inter- 
actions between dust particles, long-range and short-range inter- 
actions. The long-range interaction can be seen as a collective 
process involving interactions between clouds of particles. On 
the other hand, the short-range interaction is important for re- 
solving close approaches between pairs of particles (i.e., gravi- 
tational scattering). For the sake of clarity, we separate these two 
processes. 

Gravitational collapse timescale tqi 

We introduce a length scale A » 6 r , where 6 r is the average 
distance between neighbouring dust particles. In this limit, 
which we call the long-range interaction, the dust particle 
distribution can be approximated by a continuous density 
p, and one can define a gravitational timescale, which is 
of the order of the free fall time of the system, defined by 
tgi = 1 / VGp where G is the gravitational constant. 

Gravitational scattering timescale tq s 

Short-range gravitational interactions can be seen as an in- 
teraction between a pair of single particles that can result in 
a scattering event. As with physical collisions, the timescale 
for a gravitational interaction depends on the cross-section, 
tgs = (o"GVp«) _1 , where <x G is the gravitational cross-section 
of the particle. Using Kepler's laws, one finds that cr c ~ 
G 2 m p /v p . The gravitational cross-section is velocity depen- 
dent, which makes it qualitatively different from physical 
(billard ball) collisions between two particles. We note that 
gravitational focusing can also lower the physical collision 
timescale. 

As mentioned earlier (Sect. [TJ, the particles might also be af- 
fected by excitation from a turbulent background state on a 
timescale defined by the turbelence itself. In this paper, we 
model this excitation in a simplified way, described in Sect. 
[3] In our case, the turbulence is scale independent and has no 
timescale associated with it. 
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2.2. The real physical system 

To identify the dominant physical processes in a protoplanetary 
disc, one has to quantify and compare the relevant timescales 
as defined above. In the following discussion and in the numer- 
ical simulations, we assume Rq — 1 AU, a gas disc thickness 
of H/Rq = 0.01, and a minimum mass solar nebula (MMSN) 
with surface gas density S = 890 g cirT 2 . One usually assumes 
a solid to gas ratio of p s /p g = 0.01. We, however, assume a solid 
to gas ratio of unity. As ment ioned above, this value is j usti- 
fied by numerical simulations dJohansen et al.ll2007l l2Q09h that 
demonstrate overdensities of the order of 10 1 - 10 2 can easily 
occur because of the interaction with a turbulent gas disc, the 
streaming instability, or vertical settling. We note that signifi- 
cantly differe nt models of the solar nebula have been proposed 
(lDeschll2007l) . However, all our results can easily be scaled to 
different scenarios. 

We simplify the system by assuming that all solid compo- 
nents of the disc are in meter-sized particles. We assume that 
these boulders have a typical velocity dispersion caused by gas 
turbulence v p ~ 0.05c, ~ 30m/s, where c s is t he local sound 
speed , in accordance with numerical results dJohansen et al.l 
l2007h . Since v p <K c s , the particles tend to sediment toward the 
midplane, forming a finite thickness dust layer due to the non- 
zero velocity dispersion. 

For meter-sized boulders, the physical cross-section is <x c ~ 
1 m 2 , whereas the gravitational scattering cross-section is cr G ~ 
10~ 21 m 2 showing clearly that gravitational scattering is irrel- 
evant to the dynamics of dust particles embedded in a disc. 
However, all other timescales are roughly equivalent. Meter- 
sized boulders are weakly coupled to the g as with a stopping 
time t s ~ T e ~ Q. 1 dWeidenschillindri977ah . The physical col- 
lision timescale is r c = 7.3 Q _I and the long-range gravitational 
interaction timescales is tq\ ~ Sl~ l . 

This physical system is expected to bec ome gravitatio nally 
unstable, according to the Toomre criterion dToomrdll964l) . The 
instability occurs when the gravitational collapse timescale tqi 
is shorter than the transit time due to random particle motions 
A/\p and the orbital timescale. Assuming that the particles can 
be modelled by an isotropic gafl with a sound speed v p , the sys- 
tem becomes unstable when 
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Q = < ficrit - l. 
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In that case, the most unstable wavelength is given by 
2tt 2 G£ 



(1) 



(2) 



In the following, we compare the physical parameters from this 
section to their counterpart in numerical simulations. We first 
summarize the super-particle approach before presenting the 
scaling method. 



super-particles undergo gravitational scattering events, which 
is unphysical since these events never occur in a real system. 
Individual particle-particle collisions are not modelled. 

There are various examples where this approach is used 
successfully. For example, smoothed particle hydrodynamics 
(SPH) uses th e super-part icle approximation to simulate many 
gas molecules dLucvHl977t) . These systems are often assumed to 
be strongly collisional to ensure thermodynamical equilibrium 
inside each super-particle. One can therefore assign collective 
properties to clouds of particles such as pressure and tempera- 
ture. Another example is the evolution of galaxies. When two 
galaxies collide, individual particles (stars) will usually not un- 
dergo gravitational scattering events or physical collisions. In 
that case, the super-particle approach models a collisionless sys- 
tem in which collective dynamics is the only important physical 
process. In both cases, the super-particle approximation is a valid 
approach for simulating the system numerically, but will break 
down as soon as the system is marginally collisional. 

As an example, we consider two clouds of particles undergo- 
ing a "collision". In the strongly collisional case, the clouds will 
slowly merge and the thermodynamic variables (e.g., tempera- 
ture) will diffuse between the clouds, the particles inside each 
cloud following a random walk trajectory due to numerous col- 
lisions. On the other hand, in the collionless regime, the two 
clouds simply do not see each other because there is no short- 
range interaction present between the particles. In a marginally 
collisional regime, some particles will collide with particles of 
the other cloud, leading to a partial thermalisation of the velocity 
distribution, but some other particles will not have any collision 
at all and will follow approximately a straight line. Evidently, the 
outcome of this event cannot be described using a super-particle 
approach. 

2.4. Scaling method 

The idea of our scaling method is that one should keep all impor- 
tant timescales in a numerical simulation as close to those of the 
real physical system as possible and model all particle collisions 
explicitly. 

The numerical system consists of N nam particles in a box 
with length H, simulating a patch of the disc at a fixed radius 
Rq. The particle mass is m num = 2,H 2 /N num (H is the box size 
and one scale height). Thus, the density in the box and the long- 
range gravitational interactions are unchanged compared to the 
real system. 

The gravitational scattering cross-section of the numerical 
system is then given by 
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(3) 



For the initial surface density and velocity dispersion used in the 
simulation (Sect. [2T2J, one finds 



2.3. Super-particle approximation 

A super-particle represents many smaller particles. The dynam- 
ics of the small particles are not calculated exactly. It is assumed 
that all particles behave similarly and in a collective manner. To 
simulate gravitational interaction between super-particles, one 
has to use a softened potential. Without this, the gravitational 
scattering cross-section of super-particles becomes too large and 



1 This would be true for a strongly collisional system, but is not for- 
mally valid in the studied regime. 
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The physical collision cross-section in the simulation is 
o"c,num = na num where a n um, is the radius of the particles in 
the simulation. We derive two constraints from the physical and 
gravitational collision cross-sections: 

1 . Same mean free path in simulation and real physical systems 
That means Af num cr cnum = Ncr c , where cr c and N are the 
geometrical cross-section and the number of particles in a 
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region of size H 1 in a real disc, respectively. In other words, 
this condition ensures that the physical collision timescale 
in the simulation is exactly the same as in the real system. 

2. Negligible gravitational scattering 

When two particles approach each other, the outcome should 
be a physical collision, which means that cr c num » cr Gnum . 
The gravitational scattering timescale remains long com- 
pared to the physical collision timescale. 

The first condition places a constraint on the particle size in the 
simulation. With N = 4.7 ■ 10 18 and a particle radius of a = 1 m 
as found in a real disc within a box of volume // 3 , using the 
parameters from above finds that 



0" c , 



No- r 



N m 



0.01 4 



:AU. 



(5) 



We note that once this condition is satisfied in the initial condi- 
tions it will be automatically satisfied at all times. In simulations, 
we typically find that the collision timescale is reduced by more 
than one order of magnitude during the gravitational collapse. 

The second condition is then satisfied by changing the num- 
ber of particles. An interesting result is that if N num > 10, the 
second condition is easily satisfied if the first one is satisfied. 
We note however that <x G num depends strongly on the velocity 
dispersion. In particular, a velocity dispersion 5 times smaller 
than the initial value (as found in some simulations) leads to 
an increase by a factor of 600 in the gravitational scattering 
cross-section. Therefore, we suggest using a large safety factor 
(A^num > 10 5 ) to ensure that the second condition is always satis- 
fied, even for significantly smaller v p . This condition also allows 
us to estimate when our approach breaks down, namely when the 
number of clumps in the system is so low (N < 10 ~ 100) that 
gravitational scattering becomes important. 

Although it would be helpful as a further simplification, it is 
not possible to perfom the above mentioned simulation in two 
dimensions. In that case, the filling factor, which is defined to 
be the ratio of the volume (or area) of all particles to the total 
volume (or area), is of the order of one for the above parameters. 
We note that in 2D an increase in particle number (and decrease 
in particle size as required by Eq.|5J does not decrease the filling 
factor if the collisional lifetime remains constant. 



The interaction term F- mt is divided into components related to 
self gravity, physical collisions between particles, or drag and 
excitation forces: 



— Fg mv + FqoI ^drag ^ti 



(7) 



We solve the resulting equations of motion with a leap frog (kick 
drift kick) time stepping scheme. In the following subsections, 
we describe the numerical methods used to compute each of 
these terms and their physical relevance. 

The box size of 0.01 AU was chosen such that there are 
always several unstable modes in the box, when the system is 
pushed into the regime of gravitational collapse, as estimated by 
Eq.0 



ghost boxes 




Fig. 1. Shearing box, simulating a small patch of the protoplan- 
etary disc. The grey box represents the shearing box of inter- 
est, the dashed boxes surrounding the central box represent one 
ghost ring. 



3. Methods 

Two different kinds of simulation are considered in this paper: 

1 . Particles are assumed to be point masses and have no physi- 
cal size, the gravitational field of the particles being approx- 
imated with a smoothing length to avoid numerical diver- 
gences (see Sect. 13. It . 

2. Particles have a physical size and, therefore, no gravitational 
smoothing is required but physical collisions must be in- 
cluded. 

We refer to the particles as super-particles and scaled particles, 
respectively. 

We perform our simulations in a cubic box with shear peri- 
odic boundary conditions in the radial (x) and perdiodic bound- 
ary conditions in the azimuthal (y) and vertical ( z) directions, 
as illustrated in Fig[T] dWisdom & Tremaind[T988h . In the local 
approximation, the force per mass on each particle is a sum of 
the contributions from Hill's equations and the interaction terms. 
Hill's equation can be written as 

F m = -2Q e ; x \ p + 3£2 2 x e A - Q 2 z e z . (6) 



3.1. Self gravity 

In the A^-body problem that we consider, we have to solve 
Newton's equations of universal gravitation for a large number 
of particles N. The gravitational force on the i-th particle is given 
by 

G- ^ey, (8) 

j*i \ r ij + b) 

where b is the smoothing length used to avoid divergences in 
numerical simulations (in simulations that include physical colli- 
sions b — 0) and ey is the unit vector in the direction of the grav- 
itational force between the j-th and j-th particle. Calculating 
the gravity for each particle from each other particle results in 
O (n 2 ^ operations. To reduce the number of operations, one can 
use different approximatio ns to Eq. [8] In this p aper, we use a 
Barnes-Hut (BH) tree code dBarnes & HutHl986l) . 

The BH tree in three dimensions divides the original box into 
eight smaller cells with half the length of the original box. This 
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process is continued recursively until there is only one particle 
per cell left. The depth of the tree in a homogeneous medium is 
approximately log g N. One then calculates the total mass and the 
centre of mass for every cell at every level of the tree. 

To calculate the force acting on a particle, one starts at the 
top of the tree and descends into the tree as far as necessary to 
achieve a given accuracy. If the current cell is far away from the 
particle for which the force is calculated, then the detailed den- 
sity structure within this cell is not important. All that matters is 
the box's monopole moment (total mass and centre of mass). 
One therefore does not have to descend into this tree branch 
any further. The BH tree reduces the number of calculations to 
0(N log N). 

We use 8 rings of ghost boxes in the radial and azimuthal 
direction (FigQ]shows one ring). A ghost box is simply a shifted 
copy of the main box. The gravity on each particle is then cal- 
culated by summing over contributions from each (ghost) box. 
This setup approximates a medium of infinite horizontal extent 
and avoids large force discontinuities at the boundaries. We do 
not need any ghost boxes in the vertical direction because the 
disc is stratified. 

A finite number of ghost rings can only act as an approxima- 
tion of an infinite medium and the gravitational force will tend to 
concentrate particles horizontally in the centre of the box. Due 
to the 8 ghost rings, the asymmetry of the gravitational force be- 
tween the centre and the faces of the box is reduced by a factor 
of 20 compared to using no ghost rings at all. We note, however, 
that a small asymmetry remains in our simulations and leads 
to high er concentrat i ons in the box centre. In other simulations 
such as iTanga et alj d2004l) . this effect is not seen because the 
system is initially gravitationally unstable and integrated only 
for approximately one orbit. The system that we are interested 
in is marginal gravitationally unstable and this slight asymmetry 
will become important after several orbits. In the beginning of 
our simulations, we integrate for many orbits in order to reach a 
stable equilibrium. We then push the system into a gravitation- 
ally unstable regime (see Sect. l3.5l for details). During the stable 
phase, horizontal over-concentrations are to be expected and are 
indeed observed because of this asymmetry. However, since in a 
real protoplanetary disc the gravitational instability is considered 
to appear locally (e.g., inside a vortex or a similar structure), this 
effect of having a preferred location for gravitational instability 
is not unphysical and does not affect any conclusions made in 
this paper. We tested this by applying a linear cut-off to the grav- 
itational force at a distance of one box length (A. Toomre, private 
communication). In this case, there was no preferred location in 
the box and the simulation evolved in exactly the same way. 

3.2. Physical collisions 

Physical collisions between particles are treated in the follow- 
ing way. After each timestep, we check whether any two par- 
ticles are overlapping. Using the already existing tree structure 
from the gravity calculation, one can again reduce the compu- 
tational costs from o(n 2> ) to 0(NlogN). In these simulations, 
a collision between particles is defined as an overlap between 
two particles that are approaching each other. When a collision 
is detected, it is resolved assuming energy and momentum con- 
servation (perfectly elastic collisions). We note that the timestep 
has to be small enough such that no collision is missed and the 
overlap is always small compared to the particle size. 

The particle radius is given by Eq.|5] In order to keep the col- 
lision timescale t c close to unity, which is numerically the worst 



case scenario because all timescales are equivalent, we increase 
the radius by a factor of 4 in all simulations. 

3.3. Drag force 

Each particle feels a drag force. The background velocity of the 
gas \ g is assumed to be a steady Keplerian profile 



(9) 



Although accretion flows are turbulent, we choose to use a sim- 
ple velocity profile so that the behaviour of self-gravitating par- 
ticles can be understood in conditions that can be easily con- 
trolled. 



3.4. Random excitation 

The system of particles described above is dynamically unsta- 
ble even without self-gravity, as the coefficient of restitution and 
the stopping time do not depend on the particle velocitiefl The 
particles can either settle down in the midplane and create a 
razor-thin disc if the stopping time is too short, or they can ex- 
pand vertically forever if the stopping time is above some critical 
value and the excitation mechanism is provided only by colli- 
sions. Dust particles in accretion discs are found in the settling 
regime. However, a complete settling never occurs in accretion 
disks because the backgroun d flow is always turbu lent due to the 
Kelvin-Helmoltz instability dJohansen et al.ll2006|). the MRI. or 

r\th**r huHrnH\m^mi^ in ctiihi li ti {Qf*f* f* tr IT ^cnr Mr P^r\^lr\v 



other hydrodynamic instabilities (see e .g.. iLesur & P apaloizou 
2009 ) which diffuse particles vertically (iFromang & Papaloizou 
20061) . This stirring process of turbulence is not present in the 
gas velocity field of the simulations presented here (see Eq. [9j. 
To approximate the turbulent mixing, we added a random ex- 
citation (white noise in space and time) to the particles in the 
simulation. This allows us to have a well defined equilibrium 
in which the system is stable rather than starting from unstable 
initial conditions that might influence the final state. 

We perturb the velocity components of each particle after 
each timestep on a scale Av, = Vsi^, where 5t is the current 
timestep and £ is a random variable with a normal distribution 
around and variance s. This excitation mechanism heats up the 
particles and allows us to have a well defined three dimensional 
equilibrium as shown in Sect. [4] In our simulations, we use a 
value of s = 1.3 ■ 10 15 mV 3 and s = 8.9 ■ 10 14 mV 3 for the 
simulations without and with collisions, respectively. These val- 
ues were chosen such that the equilibrium state is approximately 
equivalent for both types of simulations. 



3.5. Initial conditions 

All particles, which have equal mass, smoothing length, and 
physical size, are placed randomly inside the box in the x - 
y plane. We note that particles have either a physical size or a 
smoothing length associated with them, depending on the type 
of simulation we perform (Sect. [3}. In the z-direction, the parti- 
cles are placed in a layer with an initially Gaussian distribution 
about the midplane and a standard deviation of 0.05 H. We allow 
the system to reach equilibrium by integrating it until t = 30Q~' 
(4.8yrs). 



2 This instability is qualitative ly similar to the instability described 
by Goldreich & Tremaine ( 1978) for Saturn rings, although in the latter 
case the drag forces are absent. 
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Fig. 2. Super-particles: Velocity dispersion as a function of time 
in four different runs. The simulation labeled LS has a ten times 
larger smoothing length. The curves do not overlap because the 
simulations are not converged. 
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Fig. 3. Scaled-particles: Velocity dispersion as a function of time 
in three different runs that include physical collisions. All curves 
overlap because the simulations are converged. 



length acts as an effective smoothing length. In general, to check 
the numerical resolution, the particle number N nnm is increased 
and the smoothing length b is reduced independently. A simu- 
lation is resolved when the result is independent of both N num 
and b. This turns out to be impossible in the present situation. 
The main reason is that the smallest scale in a gravitational 
collapse will ultimately depend on the smoothing length. By 
varying both parameters at the same time, an empirical scal- 
ing of b ~ 1 / V^Vnum works fairly well if N nnm is large enough. 
However, this procedure is not justified and is unphysical. The 
smoothing length was introduced to avoid divergent terms and 
not to model any small-scale physical process. Therefore, it 
should not have any impact on the physical outcome of the sim- 
ulation. Incidentally, the existence of this dependency indicates 
that short-range interactions are important for the result of these 
simulations and should be modelled with care. 

In Fig|2] we plot the velocity dispersion as a function of time. 
The simulations begin from the stable equilibrium described in 
Sect. [3] After t = 30 Q -1 , we switch off the excitation mecha- 
nism. Because the system continues to cool, it becomes gravita- 
tionally unstable within one orbit, as estimated by Eq. Q] Once 
clumping occurs, the velocity dispersion begins to rise again. All 
simulations use the particle radius given by Eq.[3]as a smoothing 
length except the ones labeled LS, which use a ten times larger 
value. The larger softening length is approximately a 128-th of 
the box length and illustrates the kind of evolution expected from 
an FFT method using a 128 3 grid. 

Snapshots of the particle distribution of two simulations are 
shown in Figj4] Both simulations use 160 000 particles and all 
parameters, except the smoothing length, are the same. The top 
row is the simulation with a smoothing length given by Eq. [5] 
whereas the bottom row uses a smoothing length that is ten 
times larger. The simulations correspond to the blue (medium 
dashed) and red (solid) curve in Figj2j We show the snapshots to 
illustrate the importance of the smoothing length in simulations 
without physical collisions. One can see that the simulations dif- 
fer already before clumps form. Stripy structures appear on a 
scale given by Eq. [2] As soon as we enter the unstable regime 
at t ~ 32fl -1 (i.e., the velocity dispersion begins to rise, see 
also Fig|2]) the simulations evolve very differently. The simula- 
tion on the top row of Fig|4]forms many clumps at an early time, 
whereas the bottom row simulation forms only a few, more mas- 
sive clumps at later times. 



We tested various ways of pushing the system into the unsta- 
ble regime, either by reducing the turbulence stirring or by short- 
ening the stopping time, but did not see any qualitative difference 
as long as we start from a well defined equilibrium state. In the 
simulations presented in this paper, we switch off the turbulence 
stirring. Following this modification, the system becomes gravi- 
tationally unstable and bound clumps form within a few orbits. 

4. Results 

4.1. Super particles 

We first present simulations without physical collisions (super- 
particles) which rely on the smoothing length b to avoid any 
divergencies. Although we use a tree code and therefore a 
smoothed potential for the force calculation is assumed, the 
results are equivalent to an FFT-based method where the grid 



This can be confirmed by looking at spectra of the same two 
simulations as shown in Fig (6] These spectra were generated by 
mapping the particles onto a 128 x 128 grid in the x — y plane 
and computing the fast Fourier transform of the mapping in the 
x direction. The resulting spectra were finally averaged in the 
y direction (see iTanga e t al. 2004, for a complete description 
of the procedure). As suggested by the snapshots, the nonlin- 
ear dynamics of the gravitational instability strongly depends 
on the smoothing length used. In particular, one observes that 
smaller scales are amplified more slowly for larger smoothing 
lengths. The resulting spectra at t — 40f2~' also differ signif- 
icantly. With a small enough smoothing length, the spectrum 
looks almost flat, whereas a large smoothing length introduces 
a cutoff at k/2n 10. We also note that the smoothing length in 
the latter case is of the order of the grid size (k/2n ~ 100). The 
cutoff observed clearly demonstrates that the smoothing length 
modifies the dynamics on scales up to 10 times larger than the 
smoothing length itself. 



Hanno Rein et al.: The validity of the super-particle approximation 7 




Fig. 4. Super-particles: Snapshots of the particle distribution in the xy plane. Both simulations use 160 000 particles with different 
smoothing lengths. The simulation on the bottom uses a ten times larger smoothing length than the one on the top. The snapshots 
were taken (from left to right) at t — 0, 30, 37,40Q~'. With a large smoothing length, the outcome looks very different, the system 
is more stable, more stripy structure can be seen and clumps form later, if at all. 




Fig. 5. Scaled particles: Snapshots of the particle distribution in the xy plane. The simulations (from top to bottom) use 40 000, 
320 000, 640 000 particles with their physical size given by Eq. In all simulations, the collision timescale was kept constant. The 
snapshots were taken (from left to right) at t = 0, 30, 37, 40 Q 1 . The simulation with 40 000 particles has a large filling factor by 
the last frame, which prevents clumps from forming. The intermediate resolution simulation (middle row) and the highest resolution 
simulation (bottom row) have more and smaller particles and thus a smaller filling factor. The results (i.e., number of clumps in the 
last frame) are very similar in the intermediate and high resolution simulations. 



8 



Hanno Rein et al.: The validity of the super-particle approximation 



4.2. Scaled particles 

The velocity dispersion evolution of simulations with N = 
1 60 000, 320 000, 640 000 particles including physical collisions 
is presented in Fig|3] In these runs, no gravitational smoothing 
length is needed. Again, the simulations start from a stable equi- 
librium and after t — 30 Q _1 , we change the same parameters as 
in the non-collisional case to push the system into the gravita- 
tionally unstable regime. Snapshots of the particle distributions 
for three runs are also plotted in Fig|5] The middle and bottom 
rows of snapshots correspond to the purple (medium dashed) and 
dark blue (long dashed) lines of Fig [3] 

One can see in Figj3] as well as Fig(5] that the intermedi- 
ate and high resolution runs produce very similar results. We 
call those simulations converged, as a change in particle num- 
ber does not change the outcome. An additional, lower resolu- 
tion run (N = 40 000) is not converged because the filling fac- 
tor is too large in dense regions, preventing clumps from being 
gravitationally bound. The problem does not occur in the non- 
collisional cases because particles are allowed to overlap. Since 
the filling factor scales with particle number as 1 / y/N num , this 
issue is resolved for runs with more than a few hundred thou- 
sand particles. We note that the velocity dispersion might differ 
slightly at later stages for the converged runs because the final 
clump radius is still determined by the physical particle radius 
and therefore by the particle number. 

We show the density distribution spectrum in Fig|7]as a func- 
tion of the number of particles. The resulting spectra do not de- 
pend on the number of particles in the system, as expected. In 
comparison with Fig|6j the simulations with collisions are con- 
verged since the dynamic is not controlled by the numerical pa- 
rameter N, which is the only free parameter in the system. 

At very late times (t > 38f2~'), the spectra begin to show 
a systematic trend towards more structure on smaller scales for 
runs with larger N. This is expected and due to the filling factor, 
already mentioned above. As soon as one clump in the simu- 
lation is only determined by the size of the individual particles 
it consists out of, our approach breaks down. One can also see 
these clumps forming by looking at the velocity dispersion in 
Fig[3] The particles inside a clump begin to dominate the ve- 
locity dispersion over the background after t ~ 38£2~'. A clear 
indication of this is the spiky structure with a typical correla- 
tion time of ~ Q corresponding to different clumps interacting 
and merging with each other. One way around this issue, which 
will be considered in fu ture work, is to allow particles to merge 
dMichikoshi et al.ll2009l) . Using this accretion model, the mass 
of the clump can be used as an upper limit. 



5. Discussion and implications 

In this paper, we have shown that convergence in A^-body sim- 
ulations of planetesimal formation via gravitational instability 
can be achieved when taking into account all relevant physical 
processes. It is absolutely vital to simulate gravity, damping, ex- 
citation, and physical collisions simultaneously, as the related 
timescales are of the same order and therefore all effects are 
strongly coupled. 

A set of simulations is defined to be converged when the re- 
sults do not depend on the particle number or any other artificial 
numerical parameter such as a smoothing length. As a test case, a 
box with shear periodic boundary conditions was used, contain- 
ing hundreds of thousands of self-gravitating (super-)particles 
in a stratified equilibrium state. The particles were then pushed 
into a self gravitating regime that eventually led to gravitational 



N=160000 

t=38 



N=1 60000 LS 
t=36 




' • t= 



N=160000 LS 
V ; 1=38 



N=1 60000 LS 
N=1 60000 
t=30 



1 10 100 

Wavelength 1/L 

Fig. 6. Super-particles: Density distribution spectrum with 
160 000 particles and two different smoothing lengths that differ 
by a factor 10 at times t — 30 (red, solid curve), t — 36 (green, 
long dashed curve), and t — 38 (blue, short dashed curved). The 
spectra at similar times are not on top of each other because the 
simulations are not converged. 
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Fig. 7. Scaled particles: Density distribution spectrum with 
160 000, 320 000, and 640 000 particles (colors are equivalent to 
Fig|6]i. At each time, the spectra are converged, i.e., the spectra 
are on top of each other and independent of the particle num- 
ber, besides the noise level on very small scales. At late times 
(t > 38, blue curve), one begins to see more structure on small 
scales in simulations with larger number of particles (see text). 



collapse. Simulations with and without physical collisions were 
studied for a range of particle numbers to test convergence. 

In cases without physical collisions, convergence can not be 
achieved. However, it was possible to change multiple free sim- 
ulation parameters at the same time, namely the particle number 
N and the smoothing length b, such that in special cases the re- 
sults did not depend on the particle number. We note however 
that there is a free parameter in the simulation (i.e., the smooth- 
ing length) that effects the outcome and makes it impossible to 
find the real physical solution. 



Hanno Rein et al.: The validity of the super-particle approximation 



9 



In the protoplanetary disc, physical collisions dominate over 
gravitational scattering. In the case of planetesimal formation, 
the system is marginally collisional. Physical collisions will par- 
tially randomise the particle distribution, whereas a smoothing 
length does not because the softening length is very large com- 
pared to the gravitational cross-section. Even if the gravitational 
particle-particle scattering becomes important, it can be seen 
from Eq. [3] that the velocity dependence of the cross-section 
cr G differs fundamentally from the velocity independent cross- 
section cr c . We therefore argue that the super-particle approach 
should not be used when collisions become important. 

With collisions, the simulation outcome is both quantita- 
tively and qualitatively very different from simulations with no 
physical collisions. The initial size of the clumps is larger and the 
number of clumps is smaller. As there is no smoothing length, 
there is effectively one fewer free parameter. Changing the par- 
ticle number while keeping the collision time constant does not 
change the outcome. We therefore call these simulations con- 
verged. 

Additional tests including inelastic collisions with a normal 
coefficient of restitution of 0.25 have been performed to confirm 
that the results do not depend on the way physical collisions are 
modelled. As expected, the qualitative outcome in terms of num- 
ber of clumps and clump size is different because the physical 
properties of the system have been changed. However, once in- 
dividual collisions are resolved the results converge in exactly 
the same manner as in those simulations presented above (which 
have a coefficient of restitution of 1). 

We have also explored other gravity solvers by comparing a 
fast Fourier (FFT) gravity solver with the BH tree code (using a 
smoothing length) and confirmed our previous results. The grid 
in the FFT code acts as a smoothing length and a direct com- 
parison between the two gravity solvers gave an approximate ef- 
fective smoothing length of a quarter of the grid length. Because 
one has to solve individual particle-particle interactions, the grid 
size has to be smaller than the physical size of the particles. As 
a result, the FFT method is always slower than a tree code. In 
addition, the tree structure can be reused to search efficiently for 
collisions. 

This paper focused on the numerical requirements to study 
gravitational instability and the formation of planetesimals in 
protoplanetary discs. The initial conditions were chosen such 
that the equilibrium is a well defined starting point for the con- 
vergence study. The collision, damping, collapse and orbital 
timescales are all of the same order, which is effectively the 
worst case scenario. To provide any constraints on planetesimal 
formation itself, one would have to properly simulate the turbu- 
lent background gas dynamics which is beyond the scope of this 
paper. 

However, we were able to determine the numerical resolu- 
tion needed to resolve dust particles in a protoplanetary disc. 
Assuming that one wishes to simulate dust particles realisti- 
cally, including collisions, and that the particles are uniformly 
distributed in a box of base L 2 and height H, the size of each 
particle is determined by the collision rate in the real disc (see 
Eq.E) 



flnum = J — a, (10) 

where a and N are the size and number of particles in the volume 
LrH of the real disc. The number of particles in the simulation 
A^num has to be large enough to ensure that the filling factor is 
low until clumps form. The requirement that the filling factor is 



less than unity at t — is given by 

^(anu m ) 3 -N num <C f L 2 H, (11) 
3 

where Cf is a safety factor. Although Cf < 1 would be enough 
to resolve collisions initially, it is insufficient to simulate the col- 
lapse phase. During the collapse, the filling factor rises rapidly. 
Assuming that one wishes to resolve collisions correctly when 
the particles have contracted by one order of magnitude, one has 
to include a safety factor of Cf = 10~ 3 . If the collapse occurs 
mainly in two dimensions, as in the simulations presented here, a 
factor of Cf = 10 2 is sufficient. Furthermore, one has to ensure 
that the box contains at least a few unstable modes (see Eq. [2j. 
All this together places tight numerical constraints on numerical 
simulations of planetesimal formation via gravitational instabil- 
ity. 

In future work, we will expand the discussion to include 
more physics, namely a proper treatment and feedback of the 
background gas turbulence. This will allow us to further clar- 
ify the behaviour of planetesimals in the early stages of planet 
formation and eventually test the different formation scenarios. 
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